const char *file = "12C.den";

void primp(TGraph *g); // primp the graph
double simp(int n, const double *x, const double *f); // integrate using Simpson's rule
double integr(TGraph *g); // 4*pi*int{r^2*f(r)dr}
double rms(TGraph *g); // sqrt(4*pi*int{r^4*f(r)dr})/cc

static int cc = 0; // color
void readDen(){
	ifstream f(file);
	string line;
	for(int i = 1; i--;) getline(f, line);

	TCanvas *c = new TCanvas("c", "c", 800, 600);
	TLegend *lg = new TLegend(0.7, 0.5, 0.89, 0.89);
	lg->SetBorderSize(0);
	TGraph *gra = new TGraph(); primp(gra); lg->AddEntry(gra, "A");
	double r, ra;
	int i = 0;
	while(f >> r >> ra){
		gra->SetPoint(i, r, ra);
		i++;
	} // end while
	cout << "int_gra: " << integr(gra) << " rms: " << rms(gra) << endl;

	gra->Draw("apl");
	lg->Draw();
} // end main function

void primp(TGraph *g){
  cc++;
  if(5 == cc || 10 == cc) cc++; // skip yellow and white
  g->SetLineWidth(2); g->SetLineColor(cc); g->SetLineStyle(1);
  g->SetMarkerSize(1); g->SetMarkerStyle(6); g->SetMarkerColor(cc);
  TAxis *ax = g->GetXaxis(), *ay = g->GetYaxis(); ay->CenterTitle();
  ax->SetTitle("r [fm]"); ax->CenterTitle(); ax->CenterTitle();
  ay->SetTitleOffset(1.4); ay->SetTitle("\\rho [fm^{-3}]");
  // g->Write("", TObject::kOverwrite);
} // end primp

double simp(int n, const double *x, const double *f){
  // number of individual Simpson intervals
  if(n < 2){ cout << "Simpson: n is less than 2." << endl; getchar(); }
  int k = (n-1) / 2; // f[n-1] is dropped in the case where n is even

  double sum = f[0] + f[2*k];
  for(int i = 0; i < k; i++) sum += 4.*f[2*i+1];
  for(int i = 1; i < k; i++) sum += 2.*f[2*i];

  const double h = (x[n-1] - x[0]) / (n-1);
  if(h < 0.){ cout << "Simpson: x[n-1] is less than x[0]." << endl; getchar(); }
  sum *=  h / 3.;
  // Simpson's rule is only for odd n
  // otherwise integral of f[n-1]*h should be explicitly included by trapezoidal rule
  if(n % 2 == 0) sum += (f[n-1]+f[n-2]) * h / 2.;

  return sum;
} // end function simp

static const double FOURPI = 4.*3.14159265358979323846;
double integr(TGraph *g){
	const int n = g->GetN();
	double rr[n], yy[n];
	for(int i = 0; i < n; i++){
		g->GetPoint(i, rr[i], yy[i]);
		yy[i] *= rr[i]*rr[i];
	} // end for over i
	return simp(n, rr, yy)*FOURPI;
} // end function integr
double rms(TGraph *g){
	const int n = g->GetN();
	double rr[n], yy[n];
	for(int i = 0; i < n; i++){
		g->GetPoint(i, rr[i], yy[i]);
		yy[i] *= pow(rr[i], 4);
	} // end for over i
	return sqrt(simp(n, rr, yy)*FOURPI / integr(g));
} // end function integr
